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ABSTRACT 

We are motivated to model a heat transfer to a multiple layer regime and their optimization for 
heat energy resources. Such a problem can be modeled by a porous media with different phases 
(liquid and solid). 

The idea arose of a geothermal energy reservoir which can be used by cities, e.g. Berlin. 

While hot ground areas are covered to most high populated cites, the energy resources are 
important and a shift to use such resources are enormous. 

We design a model of the heat transport via the flow of water through the heterogeneous layer 
of the underlying earth sediments. 

We discuss a multiple layer model, based on mobile and immobile zones. 

Such numerical simulations help to economize on expensive physical experiments and obtain 
control mechanisms for the delicate heating process. 

Keywords: Multiple Layer Regime, Multiple phase model, convection-diffusion reaction equa- 
tions. 
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1 INTRODUCTION 



We motivate our research on simulating novel energy resources in geothermic. 

The heat transfer in permeable and non-permeable layers are models and we simulate the tem- 
peratures in the different layers. 

Such simulations allow to predict possible energy resources to geothermal reservoirs. 

For such processes, we present a multi phase and multi-species model, see (IGeiser 20 09). 

The solver methods are fast Runge-Kutta solvers, whereas the mobile terms are convection- 
diffusion equations and are solved with splitting semi-implicit finite volume methods and char- 
acteristic methods, (|Geiser 20061 ). 

Such a sequential treatment of the partial differential equations and ordinary differential equa- 
tions allow of saving computational time, while expensive implicit Runge-Kutta methods are 
reduced to the partial operators and fast explicit Runge-Kutta methods are for the ordinary 
operators of the multi phase model. 

With various source terms we control the required concentration at the final temperature area. 
This paper is outlined as follows. 

In Section [2l we present our mathematical model based on the multiphases. In Section [3l we 
discuss discretization and solver methods with respect to their efficiency and accuracy. The 
splitting schemes are discussed in Section |4] The numerical experiments are given in Section |5] 
In Section [6l we briefly summarize our results. 

2 MATHEMATICAL MODELING 

In the model we have included the following multiple physical processes, related to the deposi- 
tion process: 

• Flow field of the fluid: Navier-Stokes equation 

• Transport system of the species: mobile and immobile phases 
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In the following we discuss the three models separately and combine all the models into a 
multiple physical model. We assume a two-dimensional domain of the apparatus with isotropic 
flow fields, see (Gobbert and Ringhofer 1998). 



2.1 Flow field 

The conservation of momentum is given by (flow field: Navier-Stokes equation) 

d 

— v + v- Vv = -Vp, mil x [0,* (1) 
at 

v(x,i) = v (x), on n, (2) 
v(x,i) = vi(x,t), onffix [0,t], (3) 

where v is the velocity field, p the pressure, v the initial velocity field and the position vector 
x = (xi,x 2 )' G fi C H 2 ' + . Furthermore, we assume that the flow is divergence free and the 
pressure is pre-defined. 

2.2 Transport systems (multi phase equations) 

We model the heat transfer as an underlying medium in the earth layers with mobile and im- 
mobile phases. Here heat transport in the fluid with different species contain of mobile and 
immobile concentrations. For such a heterogeneous media, we applied our expertise in model- 
ing multiphase transport through a porous medium. 



3 



Multi-layer Regime of a porous media 




Figure 1: Multiple layer regime of the underlying rocks and earth layers. 



In the model, we consider both absorption and adsorption taking place simultaneously and 
with given exchange rates. Therefore we consider the effect of the gas concentrations' being 
incorporated into the porous medium. 
We extend the model to two more phases: 

• Immobile phase 

• Adsorbed phase 

In Figure [2l the mobile and immobile phases of the gas concentration are shown in the macro- 
scopic scale of the porous medium. Here the exchange rate between the mobile gas concentra- 
tion and the immobile gas concentration control the flux to the medium. 
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• mobile phase 




immobile phase 



Exchange 
(mobile immobile) 



Figure 2: Mobile and immobile phase. 



In Figure |3j the mobile and adsorbed phases of the gas concentration are shown in the macro- 
scopic scale of the porous medium. To be more detailed in the mobile and immobile phases, 
where the gas concentrations can be adsorbed or absorbed, we consider a further phase. Here 
the adsorption in the mobile and immobile phase is treated as a retardation and given by a 
permeability in such layers. 



o mobile phase 




q§°°o® — Kinetically 
O O cP® Controlled Sorption 




O OsPa • adsorbed phase 



■Co © O s orbed area 



Figure 3: Mobile-adsorbed phase and immobile-adsorbed phase. 
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The model equation for the multiple phase equations are 

^ + V • Fi = gi-Ti + T itim ) + k a {-T { + T i>ad ) 

+ Kk<PT k + Qi, in fi x [0, t], (4) 

k=k(i) 

F t = wTi - D e ®VTi, (5) 
-\i,i(pT itim + ^ Kk<t>Tk,im + Qi,im, infix [0,t], (6) 

k=k(i) 

4>d t T iM = k a (Ti - T iM ) - \4>T iM + ^ Kk<t>Tk,ad + Qi,ad, in fi x [0, t], (7) 

k=k(i) 

,im,ad + Q i,im,adi i m fix[0,t], (8) 

k=k(i) 

7;(x, t) = q i0 (x), T i)0d (x, t) = 0, T i)im (x, t) = 0, T ijimM (x : t) = 0, on fi, (9) 
Ti(x,t) = T Jjl (x,t),T vld (x,t) =0,T Mm (x,t) = 0,T iiimiad (x,t) = 0, on <9fi x [0,^0) 

where the initial value is given as T ij0 and we assume a Dirichlet boundary conditions with 
the function T^x, t) sufficiently smooth, all other initial and boundary conditions of the other 
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phases are zero. 



: effective porosity [— ], 

Tj : temperature of the ith species in the underlying rock 
Tj jim : temperature of the ith species in the immobile zones of the rock 
phase [K/m 3 ], 

T^ ad : temperature of the ith species in the adsorbed zones of the rock 
phase [K/m 3 ], 

T^im,ad '■ temperature of the ith species in the immobile adsorbed zones of the rock 
phase [K/m 3 ], 

v : velocity through the rock and porous substrate (|Rouch 2006) [cm/ h}] , 
_D e W : element- specific diffusion-dispersions tensor [m 2 /h]], 
: decay constant of the zth species 
Qi : source term of the ith species [K / {m 3 h)], 
g : exchange rate between the mobile and immobile concentration [1/h], 
k a : exchange rate between the mobile and adsorbed concentration or immobile and 
immobile adsorbed concentration (kinetic controlled sorption) [1/h], 

with i = 1, . . . , M and M denotes the number of components. 
The parameters in © are further described, see also (Gei ser 20031) . 

The four phases are treated in the full domain, such that we have a full coupling in time and 
space. 

The effective porosity is denoted by <p and declares the portion of the porosities of the aquifer 
that is filled with solid grain, and we assume a nearly solid phase. The transport term is indicated 
by the Darcy velocity v, that presents the flow direction and the absolute value of the heat flux. 
The velocity field is divergence free. The decay constant of the zth species is denoted by \. 
Thereby, k(i) denotes the indices of the other species. 



3 DISCRETIZATION AND SOLVER METHODS 



We first discretize the underlying flow and transport equations in space with finite volume meth- 
ods, while we then apply the time integration methods, e.g. Runge-Kutta schemes. 

3.1 Notation 

The time-steps for the calculation in the time-intervals are (t n , t n+l ) C (0, T) , for n = 0, 1, 

The computational cells are given as Qj C with j = 1, . . . , I. The unknown I is the number 
of the nodes. 

For the application of finite- volumes we have to construct a dual mesh for the triangulation T 
, for the domain f2. First the finite-elements for the domain Cl are given by T e , e — 1, . . . , E. 
The polygonal computational cells £lj are related to the vertexes Xj of the triangulation. 

The notation for the relation between the neighbor cells and the concerned volume of each cell 
is given in the following notation. 

Let Vj = | IX, | and the set Aj denote the neighbor-point x k to the point Xj. The boundary of the 
cell j and k is denoted as T jk . 

We define the flux over the boundary T jk as 

Vj k — / n ■ v ds . (11) 

The inflow-flux is given as v jk < 0, and the outflow-flux is v jk > 0. The antisymmetry of the 
fluxes is denoted as Vj k = —v k j. The total outflow-flux is given as 

k^out(j) 

The idea of the finite- volumes is to construct an algebraic system of equation to express the 
unknowns c" c(xj,t n ). The initial values are given by c°. The expression of the interpolation 
schemes can be given naturally in two ways: the first possibility is given with the primary mesh 
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of the finite-elements 

i 

c » = ^ c -0 J (x) (13) 

where are the standard globally-finite element basis functions (Frolkovic and Geiser 2003). 
The second possibility is given with the dual mesh of the finite volumes with, 

5 n = X>>i(*) (14) 

3=1 

where cpj are piecewise constant discontinuous functions defined by (fj (x) = 1 for x G fi,- and 
<Pj(x) = otherwise. 

3.2 Discretization of the Transport equation 

We deal with the transport part, see ©: 



Ri— q + VF % = 0, in ft x [0,t] (15) 

= vq - D e «Vq, 

Cj(x,t) = Cj i0 (x), on Q, (16) 

Ci(x,t) = Cj 5 i(x, t), on <9f2 x [0,i], (17) 



For the convection part, we use a piecewise constant finite volume method with upwind dis- 
cretization, see (IFrolkovic and Geiser 20031 ). For the diffusion-dispersion part, we also apply a 
finite volume method and we assume the boundary values are denoted by n-D e ^ Vcj(x, t) = 0, 
where x 6 T is the boundary T = dVl, cf. (Frol kovic 2002al) . The initial conditions are given 

by ct(x,0) = Ci t0 (x). 

We integrate (fT5l) over space and obtain 

/ dx= ! V ■ (— vq + 7J e(i) Vci) dx . (18) 
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The time integration is done later in the decomposition method with implicit-explicit Runge- 
Kutta methods. Further the diffusion-dispersion term is lumped, cf. (|Geiser 20031) Eq. (TT8T ) is 
discretized over space using Green's formula. 

VjRi^Ci dx = ^ n • (-vq + D e «Vcj) d 1 , (19) 

where Tj is the boundary of the finite volume cell Qj and V u j is the volume of the cell j. We 
use the approximation in space, see (IGeiser 20031) . 



The spatial integration for the diffusion part (1191 ) is done by the mid-point rule over its finite 
boundaries and the convection part is done with a flux limiter and we obtain: 

V}/^|-o, = Yl n e Vvc^ 7 + E E l r i*l n i* • D h^ c hk , (20) 



where |T^ fc | is the length of the boundary element r%. The gradients are calculated with the 
piecewise finite-element function fa. 

We decide to discretize the ux with an up-winding scheme and obtain the following discretiza- 
tion for the convection part: 



Vj e Cij if V j e > 0, 
Fj, e = <( (21) 

Vj, e Ci,k if Uj, e < 0, 



where fj >e = f v • \ij^ds. 

We obtain for the diffusion part: 



V Cj e fc = E c < V<Mx, e fc ) . (22) 



We get, using difference notation for the neighbor points j and I, cf. (Frolkovic and De Schepper 2001 1, 
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the full semi-discretization: 



v;/'V-,, = E F ^ + E E ( E m n h • ^u^S) (<* - <*) . 

eeAj eGAj /GA e \{j} fcGAj 

where j = 1, . . . , m. 

Remark 1 For higher order discretization of the convection equation, we apply a reconstruc- 
tion which is based on Godunov's method. We apply a limiter function that fulfills the local 
min-max property. The method is explained in HFrolkovic and Geiser 2003\) . The linear poly- 
nomials are reconstructed by the element-wise gradient and are given by 

u(xj) = Cj , (23) 
1 E f 

Vu\v 3 = tfE / Vccte ' (24) 
with j = 1, . . . , I . 

The piecewise linear functions are denoted by 

Ujk = C j + ^jVu\ Vj (x jk - Xj) , (25) 
with j — 1, ... ,1 , 

where ipj G (0, 1) is the limiter function and based on this, d25l) fulfills the discrete minimum 
maximum property, as described in AFrolkovic and Geiser 2003\ . 

3.3 Discretization of the source-terms 

The source terms are part of the convection-diffusion equations and are given as follows: 

d tCi (x, t) - v • Vcj + VDVa = ft (x, t) , (26) 

where i = 1, . . . , m, v is the velocity, D is the diffusion tensor and ft(x,t) are the source 
functions, which can be point wise, linear in the domain. 
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The point wise sources are given as : 



qi(t) 



^ t<T, 



with / qi(t)dt = q St 



(27) 



t>T, 



r 



where q s i is the concentration of species i at source point x sourceii 6 over the whole time- 
interval. 

The line and area sources are given as : 



where q Sji is the source concentration of species i at the line or area of the source over the whole 
time-interval. 

For the finite- volume discretization we have to compute : 



where T source y is the boundary of the finite-volume cell Q sourC e,i,j which is a source area. We 
have UjQsource^j = Q SO urce,i where j e I source, where I source is the set of the finite- volume cells 
that includes the area of the source. 

The right-hand side of (l29l) is also called the flux of the sources (Frol kovic 2002bl) . 

3.4 Discretization of the Navier-Stokes equation 

We deal with the following Navier-Stokes equation: 




(28) 




(29) 



— v + v • Vv 

ut 

V-v = 0, 



Vp, in Q x [0, t] 



(30) 



(3D 
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where v = {v^v^f, for simplicity we have normalized with p — 1, and p is the pressure which 
is predefined. 



For the time discretization, we use the explicit Euler method given by: 



v 



.rt+l _ n 



Atv n • Vv" - 



AtVp n , in n 



(32) 



V-v n = 



(33) 



where At is the local time step. 

For the spatial discretization, we apply finite volume methods on staggered grids and discretize 
in each direction of the 2D Cartesian grid. The convection term in the v x-momentum equation 
is given by, see (?): 



where V h is the control volume with grid size h and S h is the underlying boundary. We integrate 
over each face of the finite volume respecting the direction of the normal vector, see (?) and 
next subsection. 

The same procedure is also used for the convection term in the v 2 momentum equation. 
3.5 Time discretization methods 

We deal with higher order time-discretization methods. We apply the Runge-Kutta methods as 
time-discretization methods to reach higher order results . 

Based on the spatial discretized transport or flow equations we obtain the following equations: 



where A is the stiffness operator and B is the reaction operator for the transport equations. f(t) 
is the right hand side, e.g. source term of the equations. 




(34) 



d t c(t) = Ac(t) + Bc(t) + f(t), < t < T , 



(35) 



c(0) = c , 
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For such a system of ordinary differential equations, we apply the Runge-Kutta methods. 

Runge-Kutta method 

We use the implicit trapezoidal rule: 







1 



i 

2 



1 
2 



(36) 



i 

2 



1 
2 



Remark 2 We apply also higher order Runge-Kutta schemes. Based on the spatial discretisa- 
tion method, which is second order finite volume schemes, we obtain the best results with second 
order RK schemes. 

4 SPLITTING METHODS 

In the following, we discuss splitting methods to decouple the system of differential equations 
to simpler parts and accelerate the solver process. 
We concentrate on two ideas: 

• Additive Splitting schemes , 

• Iterative Splitting schemes . 

4.1 Additive Splitting schemes 

We deal with the following equation: 



v P 



0=1 0=1 

u a (0) = u afi , a = 1,2, ... , p. 



a = 1,2,..., p, 



(37) 



(38) 



Further we assume A and B are self-adjoint. 
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We apply the discretization with the schemes of weights and obtain: 



„.n+l _ _.n 

B A(au n+l + (1 - a)u n ) = n , (39) 

T 

(f) n = f(at n+1 + (l-a)t n ), (40) 



By the transition to a new time level, we require: 

(B - Aar)u n+1 = <f) n , (41) 

the original problem can be transferred to 

v 

- A aP ar)u n p +1 = C a = 1, 2, . . . ,p. (42) 

By the conduction to a sequence of simpler problems we 

(B aa - ^A aa aT)u n p +1/2 = a = 1, 2, . . . ,p, (43) 
(B aa - ^A aa ar)u} +1 = a = 1, 2, . . . ,p, (44) 

(45) 

Here we have the benefit to invert only the diagonal parts of the matrices and use the idea to 
solve the triangular splitting of the operator A = A\ + A 2 . 

Theorem 1 //we choose a > \, then the splitting scheme 091) is absolute stable in an appro- 
priate Hilbert space. 

Proof 1 The outline of the proof is given in AVabishchevich 201 1\> . 
4.2 Iterative splitting method 

The following algorithm is based on the iteration with fixed- splitting discretization step-size r, 

namely, on the time-interval [t n 1 t n+l ] we solve the following sub-problems consecutively for 
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i = 0, 2, . . . 2m. (cf. (IGlowinski 20031 



Kanney and Kelley 2003).): 



dr(t) 

-^i = A lCt (t) + A 2 c l . 1 (t), with a{t n ) = c n (46) 
at 

and c (t n ) = c n , c_! = 0.0, 

^^=Aic i (t) + A 2 c m (t), (47) 
with c i+1 (t n )=c n , 



where c n is the known split approximation at the time-level t = t n . The split approximation 
at the time-level t = t n+1 is defined as c n+1 = c 2m+ i(t n+1 ). (Clearly, the function c i+1 (t) 
depends on the interval [t n , t n+1 ], too, but, for the sake of simplicity, in our notation we omit 
the dependence on n.) 

In the following we will analyze the convergence and the rate of convergence of the method 
(|46|)-(|47T) for m tends to infinity for the linear operators Ai,A 2 : X — > X, where we assume that 
these operators and their sum are generators of the C semi-groups. We emphasize that these 
operators are not necessarily bounded, so the convergence is examined in a general Banach 
space setting. 

The novelty of the convergence results are the reformulation in integral-notation. Based on 
this, we can assume to have bounded integral operators which can be estimated and given in a 
recursive form. Such formulations are known in the work of (Hansen and Ostermann 20091) and 
( Jah nke and Lubich 20091) and estimations of the kernel part with the exponential operators are 
sufficient to estimate the recursive formulations. 

4.3 Splitting Method to couple mobile and immobile and adsorbed parts 

The motivation of the splitting method are based on the following observations: 



• The mobile phase is semidiscretised with fast finite volume methods and can be stored 
into a stiffness-matrix. We achieve large time steps, if we consider implicit Runge-Kutta 
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methods of lower order (e.g. implicit Euler) as a time discretization method. 

• The immobile, adsorbed and immobile-adsorbed phases are purely ordinary differential 
equations and the each cheap to solve with explicit Runge-Kutta schemes. 

• The ODEs can be seen as perturbations and can be solved all explicit in a fast iterative 
scheme. 

For the full equation we consider the following matrix notation: 



d t c = A t c + A 2 c + Bi(c - c im ) + B 2 (c- c ad ) + Q , 

dt^-im -^2^im B\{Ci m c) ~\- B 2 {Ci m Cim,ad) Qim j 

d t c ad = A 2 c ad + B 2 (c ad - c) + Q ad , 



(48) 
(49) 
(50) 
(51) 



where c = (ci, . . . , c m ) T is the spatial discretised concentration in the mobile phase, see equa- 
tion ©, c im = (cijm, . . . , c m> j m ) T is the concentration in the immobile phase, the some also 
for the other phase concentrations. A\ is the stiffness matrix of equation ©, A 2 is the reaction 
matrix of the right hand side of ©, Bi and B 2 are diagonal matrices with the exchange of the 
immobile and kinetic parameters, see equation © and dE). 

Further Q, . . . , Q im , ad are the spatial discretised sources vectors. 

Now we have the following ordinary differential equation: 



/ 



d t C 



A 1 + A 2 + B 1 + B 2 -B l -B 2 
-B 1 A 2 + B 1 + B 2 

-B 2 A 2 + B 2 





-B 2 




\ 







-Bo 







C + Q, (52) 



A 2 + B 2 J 



where C = (c, c im , c ad , c imM ) T and the right hand side is given as Q = (Q, Q im , Q ad , Q im> 



ad) 
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For such an equation we apply the decomposition of the matrices: 



d t C = AC + Q, 

d t C = 1 : C + A 2 C + Q, 



(53) 
(54) 



where 



t A 1 + A 2 ^ 



A 2 
A 2 
A 2 



£i + 5 2 -fix -B 2 

-B 1 B l + B 2 -5 2 

-B 2 5 2 

-B 2 fi 2 



(55) 



The equation system is numerically solved by an iterative scheme: 

Algorithm 1 We divide our time interval [0, T] into sub-intervals [t n , t n+1 ], where n = 0, 1, . . . N, 

t° = and t N = T. 

We start with n = 0: 

1. ) The initial conditions are given with Co(t n+1 ) = C(t n ). We start with k = 0. 

2. ) Compute the fix point iteration scheme given as: 



d t C k = AiC k + A 2 C ft-1 + Q 



(56) 



where k is the iteration index, see ^Farago 1 2005 ). For the time integration, we apply Runge- 
Kutta methods as ODE solvers, see AHairer and Wanner 1992\) and AHairer and Wanner 1 996). 

3.) The stop criterion for the time interval [t n , t n+1 ] is given as: 



| C *( t n+l) _ C k-l( t n+lN 



< err, 



(57) 



where \ \ ■ \ \ is the maximum norm over all components of the solution vector, err is a given error 
bound, e.g. err = 10~ 4 . 
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If equation ( 1571) is fulfilled, we have the result 

C(t" +1 ) = C fc (t n+1 ), (58) 

Ifn — N then we stop and are done. 

If equation ( 1571) is not fulfilled, we do k = k + 1 and go-to 2.). 

The error analysis of the schemes are given in the following Theorem: 

Theorem 2 Let A,Be £(X) be given linear bounded operators in a Banach space £(X). We 
consider the abstract Cauchy problem: 

d t C{t) = AC{t)+BC{t), t n <t<t n+1 , (59) 
C(t n ) = C n ,forn=l,...,N, (60) 

where t\ = and the final time ist^ = T G 1R + . Then problem (l59l) a unique solution. For 
a finite steps with time size r n = t n+1 — t n , the iteration (l56l) for 
k = 1,2, ... ,q is consistent with an order of consistency G(t%). 

Proof 2 The outline of the proof is given in AGeiser 2009% 

5 NUMERICAL EXPERIMENTS 

In the following, we present to heat-flow problems. 

5.1 Two phase example 

The next example is a simplified real-life problem for a multiphase transport-reaction equation. 
We deal with mobile and immobile pores in the porous media, such simulations are given for 
heat transfers in earth layers. 

We concentrate on the computational benefits of a fast computation of the iterative scheme, 
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given with matrix exponential. 



The equation is given as: 



d t d + V • Fci = g(-ci + ci )im ) - A1C1, in Q x [0, i], (61) 

d t c 2 + V • Fc 2 = #(-c 2 + c 2 , im ) + A1C1 - A 2 c 2 , in SI x [0, f], (62) 

F = v - DV, (63) 

<9 t ci )im = g(ci - c 1>im ) - Aici )iro , in Q x [0,t], (64) 

<9 t c 2)im = 5f(c 2 - c 2 ,im) + AiCi )im - A 2 c 2jim , in Q x [0, t], (65) 

ci(x, t) = ci >0 (x), c 2 (x, t) = c 2)0 (x), on Q, (66) 

ci(x, t) = ci 5 i(x, t),c 2 (x, t) = c 2 ,i(x, t), on (9fi x [0,t], (67) 

ci )im (x, t) = 0, C2 )im (x, i) = 0, on Q, (68) 

ci )im (x, t) = 0, C2 )im (x, t) = 0, on 9fi x [0, t], (69) 



In the following we deal with the semidiscretized equation given with the matrices: 



ac 



V 



A — A 1 — G G 

Ai A - A 2 - G 

G -Ai - G 

G Ai 




G 


-A 2 -G 



C, 



(70) 



where C = (ci, c 2 , c lim , c 2im ) T , while ci = {c\ : i, . . . , cij) is the solution of the first heat 
species in the mobile phase in each spatial discretization point (i = 1, . . . , I), the same is also 
for the other solution vectors. 
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We have the following two operators for the splitting method: 



f-2 1 

1 -2 1 



A 



D 



\ 

( 1 

-1 1 



+ 



Ax 



1 -2 1 
1 -2 

\ 



where I is the number of spatial points. 



Ai 



1 1 
-1 1 



e JR 



Ixl 



( Ai 
Ai 



Ai 
Xi J 



e JR 



Ixl 



( A 2 

A 2 



A 2 = 



A 2 
A 2 j 



G JR 



Ixl 
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/ 



G 



9 

g 



V 



g 

o g J 



e IR 



(75) 



We decouple into the following matrices: 

/ 

A x = 



A 














.4 

































G IR 



47x47 



(76) 



A, 



A, 



-Ai 











Ai 


-A 2 














-Ax 











Ai 


-A 


-G 





G 








-G 





G 


G 





-G 








G 





-G 



\ 



G IR 



47x47 



G IR 



47x47 



(77) 



(78) 



For the operator Ai and A 2 = A 2 + A 3 we apply the iterative splitting method. 

Based on the decomposition, operator A\ is only tridiagonal and operator A 2 is block diagonal. 
Such matrix structure reduce the computation of the exponential operators. 

The Figure |4] present the numerical errors between the exact and the numerical solution. Here 
we obtain optimal results for one-side iterative schemes on operator B, means we iterate with 

22 



respect to B and use A as right hand side. 



Remark 3 For all iterative schemes, we can reach faster results as for the The iterative schemes 
with fast computations of the exponential matrices standard schemes. With 4 — 5 iterative steps 
we obtain more accurate results as we did for the expensive standard schemes. With one-side 
iterative schemes we reach the best convergence results. 



In the following, we present a multi-layer model in the underlying rock and assume multiple 
heat sources. The aim is to see a distribution of the heat in the upper-lying earth-layers. 



5.2 Parameters of the model equations 

In the following all parameters of the model equations (H)-® are given in Table |2] 



density 
mobile porosity 
immobile porosity 
Diffusion 
longitudinal Dispersion 
transversal Dispersion 
Retardation factor 
Velocity field 
Decay rate of the 1 st heat source 
Decay rate of the 2nd heat source 
Decay rate of the 3rd heat source 


p = 1.0 
(j) = 0.333 

0.333 
D = 0.0 
a L = 0.0 
a T = 0.00 
R = 10. 0e - 4 (Henry rate). 
v = (0.0,4.0 10~ 3 )*. 

\ AB = i io- 68 . 

Aab = 2 10~ 3 , Xbnn = 1 10~ 68 . 
\ AB = 0.25 IO" 3 , \ CB = 0.5 10~ 3 . 


Geometry (2d domain) 
Boundary 


n = [o,ioo] x [o,ioo]. 

Neumann boundary at 
top, left and right boundaries. 
Outflow boundary 
at the bottom boundary 



Table 1: Model-Parameters. 



The discretization and solver method are given as: 

For the spatial discretization method, we apply Finite volume methods of 2nd order, with the 
following parameters in Table |2l 

For the time discretization method, we apply Crank-Nicolson method (2nd order), with the 
following parameters in Table |3] 
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spatial step size 


Ax m { n 1.56, Ax ma x 2.21 


refined levels 


6 


Limiter 


Slope limiter 


Test functions 


linear test function 




reconstructed with neighbor gradients 



Table 2: Spatial discretization parameters. 



Initial time-step At init = 5 10 2 

controlled time-step At max = 1.298 10 2 , At min = 1.158 10 2 
Number of time-steps 100, 80, 30, 25 

Time- step control time steps are controlled with 

the Courant-Number CFL max = 1 

Table 3: Time discretization parameters. 

For the discretised equations are solved with the following methods, see the description in Table 
4. 



Solver 


BiCGstab (Bi conjugate gradient method) 


Preconditioner 


geometric Multi-grid method 


Smoother 


Gauss-Seidel method as smoothers for 




the Multi-grid method 


Basic level 





Initial grid 


Uniform grid with 2 elements 


Maximum Level 


6 


Finest grid 


Uniform grid with 8192 elements 



Table 4: Solver methods and their parameters. 



For the numerical experiments, we discuss the heat flow of different heat sources in the under- 
lying multiple domain regime. 

The underlying software tool is r3t, which was developed to solve discretised partial differential 
equations. We use the tool to solve transport-reaction equations, see (Fein 2004). 

5.3 Temperatur in an underlying Rock with permeable and less permeable layers 

In the following we discuss the simulation with a porous media given in Figure |5] The velocity 
is given in vertical direction, the area of the domain is [0, 100] x [0, 80]. 
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Multi-layer Regime of a porous media 



velocity 
of the 
fluid 




High perineal 
(blue) 



Low permeat 
(red) 



Heat source (line source in the deep ground) 



Figure 5: Multiple layer regime of the underlying rocks and earth layers. 



In the following Figure |6] and [7J we present an example of the concentration of three inflow 

Sources XSourceliV Source! = (30,75), X s our ce2, V Sour ce2 = (50,75) and X Source3: VSource3 = 

(70, 75). The velocity is given perpendicular in the underlying layers. 



I I ! 




Figure 6: Three inflow sources x Source i,y S ou r cei = (30,75), X Sourceli VSourcel 

= (50, 75) and 

xsource3,y Sources = (70, 75) with perpendicular velocity and 2 time- steps (initialization). 
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Figure 7: Three inflow sources x Source i,y So urcei = (30,75), x Source2 ,y Source2 = (50,75) and 
xsource3, Usource3 = (70, 75) with perpendicular velocity and 150 time-steps (end phase). 

Remark 4 The numerical experiments can also befitted to real-life experiments. The problems 
are to achieve the correct diffusion and velocity-drift coefficients. The fare field simulations, we 
obtain that the temperature derivations are centered to the middle of the high permeable layers 
( in our case the layers with high heat conduction ). Such prognostic results are important to 
allow an overview, how the heat flow is distributed in the nearer earth-layers. 

6 CONCLUSIONS AND DISCUSSIONS 

We have presented a continuous model for the multiple phases, we assumed that the heat flow 
has a fluid behavior with exchange rates to adsorbed and immobile phases based on the different 
layers. 

From the methodology side of the numerical simulations, the contributions were to decouple 
the multiphase problem into single phase problems, where each single problem can be solved 
with more accuracy. The iterative schemes allows of coupling the simpler equations and for 
each additional iterative step, we could reduce the splitting error. Such iterative methods allow 
of accelerating the solver process of multiphase problems. 

We can see in the numerical experiments a loss of the heat transfer to impermeable layer and 
strong temperature gradients within permeable layers. 
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Figure 4: Numerical errors of the one-side Splitting scheme with A (upper figure), the one-side 
Splitting scheme with B (middle figure) and the iterative schemes with 1, . . . , 6 iterative steps 
(lower figure). 
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